home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / cdrift3.zip / TEC1.C < prev    next >
C/C++ Source or Header  |  1990-02-22  |  22KB  |  470 lines
  1. /* This program is Copyright (c) 1990 David Allen.  It may be freely
  2.    distributed as long as you leave my name and copyright notice on it.
  3.    I'd really like your comments and feedback; send e-mail to
  4.    davea@vlsi.ll.mit.edu, or send us-mail to David Allen, 10 O'Moore Ave,
  5.    Maynard, MA 01754. */
  6.  
  7. /* This is the file containing all of the important functions except for
  8.    trysplit (), which splits a continent into pieces.  Also, all of the main
  9.    arrays are declared here, even a couple that are only used by functions in
  10.    tec2.c.  The array declarations are first, followed by the sequencing
  11.    function onestep () and some miscellaneous routines including the text
  12.    output routines; initialization routines and the routines that do all
  13.    the interesting stuff are last. */
  14.  
  15.  
  16. #include "const.h"
  17. #include "var.h"
  18.  
  19. #include <stdio.h>
  20.  
  21. /* These defines are used for the PostScript output section of tecst(). */
  22. #define D ((double) 432 / ((XSIZE > YSIZE) ? XSIZE : YSIZE))
  23. #define ZMAX 128
  24. #define XX(x) (108 + ((x) * D))
  25. #define YY(y) (612 - ((y) * D))
  26. #define ZZ(z) (1 - (float) ((z > ZMAX) ? ZMAX : z) / (ZMAX))
  27.  
  28.  
  29. /* The following arrays are global and most are used by functions in both
  30.    source files.  The two main ones are m and t.  Each is set up to be two
  31.    2-d arrays, where each array is the size of the whole world.  M is the
  32.    map; elements in m are indices of plates, showing which squares are
  33.    covered by which plate.  T is the topography; elements in t are altitudes. */
  34.  
  35. char m[2][MAXX][MAXY]; unsigned char t[2][MAXX][MAXY];
  36.  
  37. /* Several arrays are used by the binary blob segmenter, segment() in tec2.c.
  38.    These include r, which is used to store fragment indices; many fragments
  39.    make up one region during a segmentation.  Kid is a lookup table; fragment
  40.    k belongs to region kid[k] after a segmentation is finished.  Karea[k]
  41.    is the area of fragment k. */
  42.  
  43. char r[MAXX][MAXY], kid[MAXFRAG]; short karea [MAXFRAG];
  44.  
  45. /* The merge routine gets information from the move routine; when the move
  46.    routine puts a square of one plate on top of another plate, that information
  47.    is recorded in the merge matrix mm. */
  48.  
  49. char mm[MAXPLATE][MAXPLATE];
  50.  
  51. /* The erosion routine needs an array to store delta information; during an
  52.    erosion, the increases or decreases in elevation are summed in e and then
  53.    applied all at once to the topography. */
  54.  
  55. char e[MAXX][MAXY];
  56.  
  57. /* Several routines need temporary storage for areas and plate identifiers. */
  58.  
  59. short tarea[MAXPLATE]; short ids[MAXPLATE];
  60.  
  61. /* The plates in use are stored in this data structure.  Dx,dy are the
  62.    values to move by THIS STEP ONLY; odx,ody are the permanent move
  63.    values; rx,ry are the remainder x and y values used by newdxy() to
  64.    determine dx,dy; age is the age of the plate, in steps; area is the
  65.    area of the plate, in squares; id is the value in the m array which
  66.    corresponds to this plate; next is a pointer to the next occupied
  67.    element of the plate array. */
  68.  
  69. struct plate p [MAXPLATE];
  70.  
  71. /* The linked list header for available plates and used plates are global,
  72.    as is the step counter.  */
  73.  
  74. short pavail, phead, step;
  75.  
  76.  
  77. onestep () {
  78.    /* This is the sequencing routine called by main once per step.
  79.    It just calls the important subfunctions in order:
  80.    - trysplit   finds a plate to break up, and computes new velocities
  81.    - newdxy     computes the deltas to move each plate this step
  82.    - move       moves the plates
  83.    - merge      determines results when plates rub together
  84.    - erode      erodes the terrain, adding or subtracting at the margins
  85.    - draw       draw the resulting array once every DRAWEVERY steps
  86.    The m and t arrays are double-buffered in the sense that operations go
  87.    from m[0] to m[1] or vice-versa; src and dest determine which is which. */
  88.  
  89.    short src, dest;
  90.  
  91.    src = step % 2; dest = 1 - src;
  92.    if (rnd (100) < RIFTPCT) trysplit (src);
  93.    newdxy ();
  94.    move (src, dest);
  95.    merge (dest);
  96.    if (DOERODE) erode (dest);
  97.    if (step && !(step % DRAWEVERY)) draw (dest); }
  98.  
  99.  
  100. palloc () {
  101.    /* Allocate a plate from the array and return its index.  All the fields
  102.    of the plate are initialized to 0, except `next'.  That field is used to
  103.    link together the plate structures in use.  */
  104.  
  105.    short x;
  106.  
  107.    if (!pavail) panic ("No more objects");
  108.    x = pavail; pavail = p[x].next;
  109.    p[x].next = phead; phead = x;
  110.    p[x].area = 0; p[x].age = 0;
  111.    p[x].rx = 0; p[x].ry = 0;
  112.    p[x].odx = 0; p[x].ody = 0;
  113.    p[x].dx = 0; p[x].dy = 0;
  114.    return ((int) x); }
  115.  
  116.  
  117. pfree (n) short n; {
  118.    /* Return a plate array element to the pool of available elements.
  119.       To check for infinite loops, the variable guard is incremented
  120.       at each operation; if the number of operations exceeds the maximum
  121.       possible number, the program panics. */
  122.  
  123.    short i, guard = 0;
  124.  
  125.    if (phead == n) phead = p[n].next;
  126.    else {
  127.       for (i=phead; p[i].next!=n; i=p[i].next)
  128.          if (++guard > MAXPLATE) panic ("Infinite loop in pfree");
  129.       p[i].next = p[n].next; }
  130.    p[n].next = pavail; pavail = n; }
  131.  
  132.  
  133. tecst (src, drawmode) short src, drawmode; {
  134.    /* This function is called whenever map output is called for.  It looks
  135.    at the parameter `drawmode' to decide between long text, simple text,
  136.    and PostScript output formats.  Note that the default for this
  137.    function is no output at all, corresponding to DRAWMODE_NONE. */
  138.   
  139.    register short i,j,k;
  140.  
  141.    if (drawmode == DRAWMODE_GENERIC) 
  142.       for (j=0; j<YSIZE; j++) for (i=0, k=0; i<XSIZE; i++) {
  143.          printf ("%4d", t[src][i][j]);
  144.          if (!(++k % 18)) printf ("\n"); }
  145.    else if (drawmode == DRAWMODE_TEXT) {
  146.       for (j=0; j<YSIZE; j++) {
  147.          for (i=0; i<XSIZE; i++) {
  148.             if ((k = t[src][i][j]) < ZCOAST) k = 0;
  149.             else if (k > ZMOUNTAIN) k = 2; else k = 1;
  150.             printf ("%d", k); }
  151.          printf ("\n"); }
  152.       printf ("\n"); }
  153.    else if (drawmode == DRAWMODE_GRAY) {
  154.       printf ("%%!PS-Adobe-1.0\n/d %4.2f def\n", D);
  155.       printf ("37.5 45 { dup mul exch dup mul add 1 exch sub } setscreen\n");
  156.       printf ("/r { setgray moveto d 0 rlineto 0 d rlineto ");
  157.       printf ("d neg 0 rlineto fill } bind def\n");
  158.       printf ("%6.2f %6.2f moveto\n",         XX(-1),    YY(-1));
  159.       printf ("%6.2f %6.2f lineto\n",         XX(XSIZE), YY(-1));
  160.       printf ("%6.2f %6.2f lineto\n",         XX(XSIZE), YY(YSIZE));
  161.       printf ("%6.2f %6.2f lineto\n",         XX(-1),    YY(YSIZE));
  162.       printf ("%6.2f %6.2f lineto\nstroke\n", XX(-1),    YY(-1));
  163.    
  164.       for (j=0; j<YSIZE;j++) for (i=0; i<XSIZE; i++)
  165.          if ((k = t[src][i][j]) > ZCOAST)
  166.             printf ("%6.2f %6.2f %5.4f r\n", XX(i), YY(j), ZZ(k));
  167.       printf ("\nshowpage\n"); } }
  168.  
  169.  
  170. init (s) char *s; {
  171.    /* This is the catchall function that initializes everything.  First,
  172.    it calls getparams() in tec3.c to allow the user to set parameters.  Next,
  173.    it links together the plates onto the free list and starts the used list
  174.    at empty.  The first plate is created by a fractal technique and then
  175.    improved.  Finally, the fractal is copied to the data array and drawn.
  176.    There are two kinds of improvement done here.  First, islands are
  177.    eliminated by segmenting the blob and erasing all the regions except
  178.    for the biggest.  Second, oceans inside the blob (holes) are eliminated
  179.    by segmenting the _ocean_ and filling in all regions except the biggest. */
  180.  
  181.    short besti, x; register short i, j;
  182.  
  183.    if (s) if (*s) getparams (s);
  184.  
  185.    for (i=1; i<MAXPLATE; i++) p[i].next = i + 1;
  186.    p[MAXPLATE-1].next = 0;
  187.    pavail = 1; phead = 0;
  188.  
  189.    /* Allocate a plate structure for the first plate and make a blob */
  190.    x = palloc (); makefrac (0, x);
  191.  
  192.    /* Segment m[0] looking for x, set besti to the largest region, */
  193.    /* and zero out all the other regions.  This eliminates islands. */
  194.    besti = singlefy (0, x);
  195.    if (besti > 0) for (i=1; i<XSIZE; i++) for (j=1; j<YSIZE; j++)
  196.       if (kid[r[i][j]] != besti) m[0][i][j] = 0;
  197.  
  198.    /* Segment m[0] looking for 0 (ocean), set besti to the largest region, */
  199.    /* and fill in all the other regions.  This eliminates holes in the blob. */
  200.    besti = singlefy (0, 0);
  201.    if (besti > 0) for (i=1; i<XSIZE; i++) for (j=1; j<YSIZE; j++)
  202.       if (kid[r[i][j]] != besti) m[0][i][j] = x;
  203.  
  204.    /* Fill the topo structure with the blob shape while finding its area */
  205.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++)
  206.       if (m[0][i][j]) { t[0][i][j] = ZINIT; p[x].area++; }
  207.  
  208.    /* Draw the blob */
  209.    draw (0); }
  210.  
  211.  
  212. makefrac (src, match) short src, match; {
  213.    /* This function uses a very simple fractal technique to draw a blob.  Four
  214.    one-dimensional fractals are created and stored in array x, then the
  215.    fractals are superimposed on a square array, summed and thresholded to
  216.    produce a binary blob.  Squares in the blob are set to the value in `match'.
  217.    A one-dimensional fractal of length n is computed like this.  First,
  218.    set x [n/2] to some height and set the endpoints (x[0] and x[n]) to 0.
  219.    Then do log-n iterations.  The first iteration computes 2 more values:
  220.    x[n/4] = average of x[0] and x[n/2], plus some random number, and
  221.    x[3n/4] = average of x[n/2] and x[n], plus some random number.  The second
  222.    iteration computes 4 more values (x[n/8], x[3n/8], ...) and so on.  The
  223.    random number gets smaller by a factor of two each time also.
  224.  
  225.    Anyway, you wind up with a number sequence that looks like the cross-section
  226.    of a mountain.  If you sum two fractals, one horizontal and one vertical,
  227.    you get a 3-d mountain; but it looks too symmetric.  If you sum four,
  228.    including the two 45 degree diagonals, you get much better results. */
  229.  
  230.    register short xy, dist, n, inc, i, j, k; char x[4][65];
  231.    short xofs, yofs, xmin, ymin, xmax, ymax;
  232.  
  233.    /* Compute offsets to put center of blob in center of the world, and
  234.       compute array limits to clip blob to world size */
  235.    xofs = (XSIZE - 64) >> 1; yofs = (YSIZE - 64) >> 1;
  236.    if (xofs < 0) { xmin = -xofs; xmax = 64 + xofs; }
  237.    else { xmin = 0; xmax = 64; }
  238.    if (yofs < 0) { ymin = -yofs; ymax = 64 + yofs; }
  239.    else { ymin = 0; ymax = 64; }
  240.  
  241.    for (xy=0; xy<4; xy++) {
  242.       /* Initialize loop values and fractal endpoints */
  243.       x [xy] [0] = 0; x [xy] [64] = 0; dist = 32;
  244.       x [xy] [32] = 24 + rnd (16); n = 2; inc = 16;
  245.  
  246.       /* Loop log-n times, each time halving distance and doubling iterations */
  247.       for (i=0; i<5; i++, dist>>=1, n<<=1, inc>>=1)
  248.          for (j=0, k=0; j<n; j++, k+=dist)
  249.             x[xy][k+inc] = ((x[xy][k]+x[xy][k+dist])>>1) + rnd (dist) - inc; }
  250.  
  251.    /* Superimpose fractals; if sum is greater than BLOBLEVEL, there will be */
  252.    /* land there.  x[0] is horizontal, x[1] vertical, x[2] diagonal from */
  253.    /* top left to bottom right, x[3] diagonal from TR to BL. */
  254.    for (i=xmin; i<xmax; i++) for (j=ymin; j<ymax; j++) {
  255.       k = x[0][i] + x[1][j] + x[2][(i - j + 64) >> 1] + x[3][(i + j) >> 1];
  256.       if (k > BLOBLEVEL) m[src][i+xofs][j+yofs] = match; } }
  257.  
  258.  
  259. singlefy (src, match) short src, match; {
  260.    /* This is a subfunction of init() which is called twice to improve the
  261.    fractal blob.  It calls segment() and then interprets the result.  If
  262.    only one region was found, no improvement is needed; otherwise, the
  263.    area of each region must be computed by summing the areas of all its
  264.    fragments, and the index of the largest region is returned. */
  265.  
  266.    short i, reg, frag, besti, besta;
  267.  
  268.    segment (src, match, &frag, ®);
  269.    if (reg == 1) return (-1); /* No improvement needed */
  270.  
  271.    /* Initialize the areas to zero, then sum frag areas */
  272.    for (i=1; i<=reg; i++) tarea[i] = 0;
  273.    for (i=1; i<=frag; i++) tarea [kid[i]] += karea [i];
  274.  
  275.    /* Pick largest area of all regions and return it */
  276.    for (i=1, besta=0, besti=0; i<=reg; i++)
  277.       if (besta < tarea[i]) { besti = i; besta = tarea[i]; }
  278.    return ((int) besti); }
  279.  
  280.  
  281. newdxy () {
  282.    /* For each plate, compute how many squares it should move this step.
  283.    Multiply the plate's basic movement vector odx,ody by the age modifier
  284.    MR[], then add the remainders rx,ry from the last move to get some large
  285.    integers.  Then turn the large integers into small ones by dividing by
  286.    REALSCALE and putting the remainders back into rx,ry.  This function also
  287.    increases the age of each plate, but doesn't let the age of any plate
  288.    go above MAXLIFE.  This is done to make sure that MR[] does not need to
  289.    be a long vector. */
  290.  
  291.    register short i, a;
  292.  
  293.    for (i=phead; i; i=p[i].next) {
  294.       a = (p[i].odx * MR[p[i].age]) + p[i].rx;
  295.       p[i].dx = a / REALSCALE; p[i].rx = a % REALSCALE;
  296.       a = (p[i].ody * MR[p[i].age]) + p[i].ry;
  297.       p[i].dy = a / REALSCALE; p[i].ry = a % REALSCALE;
  298.       if (p[i].age < MAXLIFE-1) (p[i].age)++; } }
  299.  
  300.  
  301. move (src, dest) short src, dest; {
  302.    /* This function moves all the plates that are drifting.  The amount to
  303.    move by is determined in newdxy().  The function simply steps through
  304.    every square in the array; if there's a plate in a square, its new location
  305.    is found and the topography is moved there.  Overlaps between plates are
  306.    detected and recorded so that merge() can resolve the collision; mountain
  307.    growing is performed.  If two land squares wind up on top of each other,
  308.    folded mountains are produced.  If a land square winds up where ocean was
  309.    previously, that square is the leading edge of a continent and grows a
  310.    mountain by subsuming the ocean basin. */
  311.  
  312.    register short i, j; short a, b, c, x, y;
  313.  
  314.    /* Clear out the merge matrix and the destination arrays */
  315.    for (i=1; i<MAXPLATE; i++) for (j=1; j<MAXPLATE; j++) mm[i][j] = 0;
  316.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
  317.       m[dest][i][j] = 0; t[dest][i][j] = 0; }
  318.  
  319.    /* Look at every square which belongs to a plate */
  320.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) if ((a = m[src][i][j]) > 0) {
  321.  
  322.       /* Add the plate's dx,dy to the position to get the square's new */
  323.       /* location; if it is off the map, throw it away */
  324.       x = p[a].dx + i; y = p[a].dy + j;
  325.       if ((x >= XSIZE) || (x < 0) || (y >= YSIZE) || (y < 0)) p[a].area--;
  326.       else { /* It IS on the map */
  327.  
  328.          /* If the destination is occupied, remove the other guy but */
  329.          /* remember that the two plates overlapped; set the new height */
  330.          /* to the larger height plus half the smaller. */
  331.          if (c = m[dest][x][y]) {
  332.             (mm[a][c])++; (p[c].area)--;
  333.             b = t[src][i][j]; c = t[dest][x][y];
  334.             t[dest][x][y] = (b > c) ? b + (c>>1) : c + (b>>1); }
  335.  
  336.          /* The destination isn't occupied.  Just copy the height. */
  337.          else t[dest][x][y] = t[src][i][j];
  338.  
  339.          /* If this square is over ocean, increase its height. */
  340.          if (t[src][x][y] < ZCOAST) t[dest][x][y] += ZSUBSUME;
  341.  
  342.          /* Plate A now owns this square */
  343.          m[dest][x][y] = a; } } }
  344.  
  345.  
  346. merge (dest) short dest; {
  347.    /* Since move has set up the merge matrix, most of the work is done.  This
  348.    function calls bump once for each pair of plates which are rubbing; note
  349.    that a and b below loop through the lower diagonal part of the matrix.
  350.    One subtle feature is that a plate can bump with several other plates in
  351.    a step; suppose that the plate is merged with the first plate it bumped.
  352.    The loop will try to bump the vanished plate with the other plates, which
  353.    would be wrong.  To avoid this, the lookup table lut is used to provide
  354.    a level of indirection.  When a plate is merged with another, its lut
  355.    entry is changed to indicate that future merges with the vanished plate
  356.    should be applied to the plate it has just been merged with. */
  357.  
  358.    char lut[MAXPLATE]; short a, aa, b, bb, i;
  359.  
  360.    for (a=1; a<MAXPLATE; a++) lut[a] = a;
  361.    for (a=2; a<MAXPLATE; a++) for (b=1; b<a; b++) if (mm[a][b] || mm[b][a]) {
  362.       aa = lut [a]; bb = lut[b];
  363.       if (aa != bb) if (bump (dest, aa, bb)) {
  364.          lut[aa] = bb;
  365.          for (i=1; i<MAXPLATE; i++) if (lut[i] == aa) lut[i] = bb; } } }
  366.  
  367.  
  368. bump (dest, a, b) short dest, a, b; {
  369.    /* Plates a and b have been moved on top of each other by some amount;
  370.    alter their movement rates for a slow collision, possibly merging them.
  371.    The collision "strength" is a ratio of the area overlap (provided by
  372.    move ()) to the total area of the plates involved.  A fraction of each
  373.    plate's current movement vector is subtracted from the movement vector
  374.    of the other plate.  If the two vectors are now within some tolerance
  375.    of each other, they are almost at rest so merge them with each other. */
  376.  
  377.    double maa, mab, ta, tb, rat, area; register short i, j, x;
  378.  
  379.    /* Find a ratio describing how strong the collision is */
  380.    x = mm[a][b] + mm[b][a]; area = p[a].area + p[b].area;
  381.    rat = x / (MAXBUMP + (area / 20)); if (rat > 1.0) rat = 1.0;
  382.  
  383.    /* Do some math to update the move vectors.  This looks complicated */
  384.    /* because a plate's actual movement vector must be multiplied by */
  385.    /* MR[age], and because I have rewritten the equations to maximize */
  386.    /* use of common factors.  Trust me, it's just inelastic collision. */
  387.    maa = p[a].area * MR[p[a].age]; mab = p[b].area * MR[p[b].age];
  388.    ta = MR[p[a].age] * area;
  389.    p[a].odx = (p[a].odx * maa + p[b].odx * mab * rat) / ta;
  390.    p[a].ody = (p[a].ody * maa + p[b].ody * mab * rat) / ta;
  391.    tb = MR[p[b].age] * area;
  392.    p[b].odx = (p[b].odx * mab + p[a].odx * maa * rat) / tb;
  393.    p[b].ody = (p[b].ody * mab + p[a].ody * maa * rat) / tb;
  394.  
  395.    /* For each axis, compute the remaining relative velocity.  If it is */
  396.    /* too large, return without merging the plates */
  397.    if (ABS (p[a].odx*MR[p[a].age] - p[b].odx*MR[p[b].age]) > BUMPTOL) return(0);
  398.    if (ABS (p[a].ody*MR[p[a].age] - p[b].ody*MR[p[b].age]) > BUMPTOL) return(0);
  399.  
  400.    /* The relative velocity is small enough, so merge the plates.  Replace */
  401.    /* all references to a with b, free a, and tell merge() a was freed. */
  402.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++)
  403.       if (m[dest][i][j] == a) m[dest][i][j] = b;
  404.    p[b].area += p[a].area; pfree (a);
  405.    return ((int) a); }
  406.  
  407.  
  408. /* The following is defined as a macro for efficiency; erosion is still
  409.    the slowest function in the simulation.  ERODE is called four times per
  410.    pixel by erode. A and b are altitudes of two adjacent squares, while c
  411.    and d are the corresponding delta values for those squares.  The amount
  412.    each square loses to an adjacent but lower square in each step is
  413.    one-eighth the difference in altitude.  This is coded as a shift right 3
  414.    bits, but since -1 >> 3 is still -1, the code must be repeated to avoid
  415.    negative deltas. */
  416. #define ERODE(a,b,c,d) \
  417. if (((a > ZSHELF) || (b > ZSHELF))) { \
  418.    if ( a > b) { x = (a - b + ERODERND) >> 3; c -= x; d += x; } \
  419.    else { x = (b - a + ERODERND) >> 3; c += x; d -= x; } }
  420.  
  421.  
  422. erode (dest) short dest; {
  423.    /* This function takes the topography in t[dest] and smooths it, lowering
  424.    mountains and raising lowlands and continental margins.  It does this by
  425.    stepping across the entire array and doing a computation once for each
  426.    pair of 8-connected pixels.  The computation is done by ERODE, above.
  427.    The computation result for a pair is a small delta for each square, which
  428.    is summed in the array e.  When the computation is finished, the delta
  429.    is applied; if this pushes an ocean square high enough, it is added to
  430.    an adjacent plate if one can be found.  Also, if a land square is eroded
  431.    low enough, it is turned into ocean and removed from its plate. */
  432.  
  433.    register short i, j, t1, x, z; short ii, jj, xx;
  434.  
  435.    /* Zero out the array for the deltas first */
  436.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) e[i][j] = 0;
  437.  
  438.    /* Step across the entire array; each pixel is adjacent to 8 others, and */
  439.    /* it turns out that if four pairs are considered for each pixel, each */
  440.    /* pair is considered exactly once.  This is important for even erosion */
  441.    for (i=1; i<XSIZE; i++) for (j=1; j<YSIZE; j++) {
  442.       t1 = t[dest][i][j];
  443.       ERODE (t1, t[dest][i][j-1],   e[i][j], e[i][j-1])
  444.       ERODE (t1, t[dest][i-1][j-1], e[i][j], e[i-1][j-1])
  445.       ERODE (t1, t[dest][i-1][j],   e[i][j], e[i-1][j])
  446.       if (j < YSIZE-1) {
  447.          ERODE (t1, t[dest][i-1][j+1], e[i][j], e[i-1][j+1]) } }
  448.  
  449.    /* Now go back across the array, applying the delta values from e[][] */
  450.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
  451.       z = t[dest][i][j] + e[i][j]; if (z < 0) z = 0; if (z > 255) z = 255;
  452.  
  453.       /* If the square just rose above shelf level, look at the four */
  454.       /* adjacent squares.  If one is a plate, add the square to that plate */
  455.       if ((z >= ZSHELF) && (t[dest][i][j] < ZSHELF)) {
  456.          if (i > 1)      if (xx = m[dest][i-1][j]) { ii = i-1; jj = j; }
  457.          if (i < XSIZE-1) if (xx = m[dest][i-1][j]) { ii = i+1; jj = j; }
  458.          if (j > 1)      if (xx = m[dest][i][j-1]) { ii = i; jj = j-1; }
  459.          if (j < YSIZE-1) if (xx = m[dest][i][j+1]) { ii = i; jj = j+1; }
  460.          p[xx].area++; m[dest][i][j] = xx; }
  461.  
  462.       /* If the square is lower than shelf level but belongs to a plate, */
  463.       /* remove it from the plate */
  464.       if (((t[dest][i][j] = z) < ZSHELF) && (x = m[dest][i][j])) {
  465.          p[x].area--; m[dest][i][j] = 0; } } }
  466.  
  467.  
  468.  
  469.  
  470.